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We obtain an exact Green-Kubo type linear response result for the heat current in an open system. The 
result is derived for classical Hamiltonian systems coupled to heat baths. Both lattice models and fluid systems 
are studied and several commonly used implementations of heat baths, stochastic as well as deterministic, are 
considered. The results are valid in arbitrary dimensions and for any system sizes. Our results are useful for 
obtaining the linear response transport properties of mesoscopic systems. Also we point out that for systems 
with anomalous heat transport, as is the case in low-dimensional systems, the use of the standard Green-Kubo 
formula is problematic and the open system formula should be used. 

PACS numbers: 



The Green-Kubo formula |lll,|2|] is a cornerstone of the study 
of transport phenomena. For a system governed by Hamil- 
tonian dynamics, the currents that flow in response to small 
applied fields can be related to the equilibrium correlation 
functions of the currents. For the case of heat transport the 
Green-Kubo formula (in the classical limit, which this paper 
is restricted to) gives: 
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K = lim lim 



dt{J{t)J{0)) , (1) 



where k is the thermal conductivity of a d-dimensional system 
of linear dimension L at temperature T. The autocorrelation 
function on the right hand side is evaluated in equilibrium, 
without a temperature gradient. The total heat current in the 
system is J{t) — /i(x, i)(ix, where j'(x, t) as the heat flux 
density. The order of the limits in Eq. ([T]l is important. With 
the correct order of limits, one can calculate the correlation 
functions with arbitrary boundary conditions and apply Eq. ^ 
to obtain the response of an open system with reservoirs at the 
ends. There have been a number of derivations of Eq. ^ by 
various authors fll Htl ■ 

There are several situations where the Green-Kubo formula 
in Eq. ([T]l is not applicable. For example, for the small struc- 
tures that are studied in mesoscopic physics, the thermody- 
namic limit is meaningless, and one is interested in the con- 
ductance of a specific finite system. Secondly, in many low di- 
mensional systems, heat transport is anomalous and the ther- 
mal conductivity diverges f?]. In such cases it is impossi- 
ble to take the limits as in Eq. ([TJ; one is there interested in 
the thermal conductance as a function of L instead of an L- 
independent thermal conductivity. The usual procedure that 
has been followed in the heat conduction literature is to put a 
cut-off at tc ^ L, in the upper limit in the Green-Kubo integral 
There is no rigorous justification of this assumption. A re- 
lated case is that of integrable systems, where the infinite time 
limit of the correlation function in Eq. ([U is non-zero. An- 
other way of using the Green-Kubo formula for finite systems 
is to include the infinite reservoirs also while applying the for- 
mula and this was done, for example, by Allen and Ford (3] 
for heat transport and by Fisher and Lee [6] for electron trans- 
port. Both these cases are for non-interacting systems and the 
final expression for conductance is what one also obtains from 



the nonequilibrium Green's function approach, a formalism of 
transport commonly used in the mesoscopic literature. More 
recently, it has been shown that Green-Kubo like expressions 
for finite open systems can be derived rigorously by using the 
steady state fluctuation theorem (SSFT)l7i|8il9[|Toll. 

In this paper, we derive a Green-Kubo like formula for open 
systems, without invoking the SSFT. Our proof applies to all 
classical systems, of arbitrary size and dimensionality, with 
a variety of commonly used implementations of heat baths. 
The proof consists in first solving the equation of motion for 
the phase space probability distribution to find the 0{AT) 
correction to the equilibrium distribution function. The aver- 
age current at this order can then be expressed in terms of the 
equilibrium correlation {J{t)Jfp{0)), where Jjp is a speci- 
fied current operator. Secondly we use the energy continu- 
ity equations to relate two different current-current correla- 
tion functions, namely ( J(0) J(t)} and ( J(0) Jf,(i)} where Jf, 
is an instantaneous current operator involving heat flux from 
the baths. Finally one relates {J{0)Jb{t)) to {J{0)Jfp(t)) 
and then, using time-reversal invariance, to ( J(i) J/p(0)}. For 
baths with stochastic dynamics, time-reversal invariance fol- 
lows from the detailed balance principle, which is an essential 
requirement of our proof. 

We first give a proof of our linear response result for a ID 
lattice model with white noise Langevin baths. We consider 
the following general Hamiltonian: 
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J2 U{xi - xi+i) , (2) 
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where x — {xj}, v = {vi] with I = \,2...N denotes dis- 
placements of the particles about their equilibrium positions 
and their velocities, and {rai} denotes their masses. The par- 
ticles at the ends are connected to two white noise heat baths 
of temperatures Tl and respectively. The equations of mo- 
tion of the system are given by: 

mvi = fi - Si,i[-/^vi - Tj^] - Si^nIi'^vn - v'^], (3) 

where // = —dH/dxi, and rj^'^{t) are Gaussian noise terms 
with zero mean and satisfying the fluctuation dissipation rela- 
tions: (7y^^^(t)77^^«(i')) = 2j^'HBTL,B5{t - t') . 
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In the first part of the proof we obtain an expression for the 
nonequilibrium steady state average (J) at, at linear order in 
AT, and then we will relate this to the equilibrium correla- 
tion function (J(i) J(0)} Corresponding to the stochas- 
tic Langevin equations in Eq. (|3]l, one has a Fokker-Planck 
(FP)equation for the phase space distribution P(x, v, t). Set- 
ting Tl^T + at/2 and Tr^T - AT/2 we write the FP 
equation in the following form; 



dt 

where T(x, v 
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where = d/dxi + {fi/mi) d/dvi ] is the 

Hamiltonian Liouville operator and 7^ = 7^,7^ 7^- 
For AT = the steady state solution of the FP equation is 
known and is just the usual equilibrium Boltzmann distribu- 
tion Pq = e^^^ /Z, where Z = J dxdve~'^^ is the canoni- 
cal partition function [(3 = (/cbT)"^]. It is easily verified that 
LPq = 0. For AT ^ 0, we solve Eq. ^ by perturbation the- 
ory, starting from the equilibrium solution at time t ~ —00. 
Writing P(x, v, t) = Po +p(x, v, <), we obtain the following 
solution at 0(AT): 
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To define the current operator, one first defines the local 
energy density at the Tth site: e; ~ mivf/2 + V{xi) + 
^[U{xi-i — xi) + U{xi — xi+i)]. Taking a time derivative 
gives the energy continuity equation 



dti/dt + ii + i,i - jij^i = Si^i + jN,R 5i^N 

where ji+ij = ^(w/ + 



(7) 



gives the current from the Tth to the I + I'th site {f 1+1,1 is the 
force on I + I'th particle due to Tth particle). We define the 

total current flowing through the system as J = X]/=i^i'+i,' ■ 
The expectation value of the total current is then given by: 



{J)/\T = j dxdv J p{x,v) 

= A/3 / dt [ dxdv J e^* Jfp Pq 



A/3 / dt{Jit)Jfp{0)) 

JQ 



(8) 



There are two parts of the proof remaining. Let us define 
the current variable as the mean of the instantaneous heat 
currents flowing into the system from the left reservoir and 
flowing out of the system to the right reservoir. Thus we have 



(9) 



where = --/^vl{t) + 7]^ {t)vi{t) , 

jNAt) = —/''v%{t)+r,''{t)vNit) . (10) 



(4) The two remaining steps then consist of proving the relations: 



(J(O)Jfc(t)) = {J{0)Jfp{t)) 

= -{j{t)jfpm, 



(11) 



(5) and / dt{J{t)J{0)) = {N - I) / dt{J{0)Jbit)) .(12) 



The first line in Eq. ( fTTT i follows from (J(0)} = and the 
result: 



{J{0)v'^{t)v,{t))^{J{0)r^"it)v^{t))=0, 



(13) 



which can be proved by making use of Novikov's theo- 
rem lfl2[ [13I1 . The second line in Eq. ( fTTT i is a state- 
ment of time-reversal symmetry. To prove this we write 
{Jfpit)Jm = JdqJdq'Jfp{q)Jiq')Po{q')Wiq,t\q',0) 
where W{q,t\q' ,0) denotes the transition probability from 
q' = (x',v') to g = (x, v) in time t. Then, using the 
detailed balance principle W{x,v,t\x'^V MlP^hc' ,v') — 
W{x', -v', t\x, -V, 0)Po(x, -v) (see H El [M]) and the 
fact that J is odd in the velocities while Jfp is even, one gets 
{Jfp{t)J{0)) = —{Jfp{0)J{t)). A more direct but equivalent 
proof is given in [17]. 

We next prove the relation given by Eq. (fT2] i. Let us define 

Diit) = ELi efc - for 1 = 1,2, ...N - 1. Then 

from the continuity equation Eq. (I?) one can show that 



dDi/dt = -2ji+i,i{t) + 2Jb{t) 



(14) 



We multiply this equation by J(0), take a steady state aver- 
age, and integrate over time from t = to 00. Since DiJ has 
an odd power of velocity we therefore get (Z?;(0)J(0)) = 
0. Also (A(oo)J(O)) = (A(oo)} (J(0)} = and 
using these we immediately get dt{ji+i,i{t)J{0)) = 
dt{Jb{t)J{0)). Summing over all bonds thus proves 
Eq. (fT2] i. With a temperature difference AT between the 
reservoirs, the steady state current between the reservoirs and 
the system {I) at is equal to (j) where j = J/{N — 1). Using 
Eqs. (1811 1112b . the conductance is given by: 
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which is the central result of the paper. 

The above proof can be extended to the case where the noise 
from the baths is exponentially correlated in time |18]. Here 
we will outline the proof for two other models: a determin- 
istic bath model coupled to a lattice Hamiltonian and another 
model where Maxwell baths are coupled to a fluid system. 
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Nose-Hoover baths: In this case the equations of motion 
are given by: mivi = fi - Si^iClVi - Si^nCrVn where (l.a 
are themselves dynamical evolving by the equations: 

Cr = iPRrriNvlf - 1)/0r ■ 

For small AT, we then write an equation of motion for the 
extended distribution function P(x, v, (j^, (fi, t) and find that 
this has the same form as Eq. (|4]l but now with: 



AT f mivl d ruNvj^ d 



2kBT^ V da Or OCr 



(16) 



If Tl = Tr = T, one can verify that the 
equilibrium phase space density is given by Pq = 
c Po(x,v) exp[— 6'lC£/2toi — dR(jj/2mN] , where c is a 
normalization constant (independent of Tl,r), and we assume 
convergence to this distribution. Acting with L^^ on this, we 
then obtain: 



Cr). 



(17) 



On the other hand, since — Cl^i is the force from the left reser- 
voir on the first particle, hence ji ^ = ^Cl^I and similarly, 
jN,R = —(rVr- Hence from the definition of in Eq. 
we obtain Jfp ^ —Ji,. The rest of the proof is similar to the 
previous case, except that there is a minus sign in the right 
hand side of the first line of Eq. ( fTTI ). This minus sign is not 
reversed in the second line of Eq. ( fTTI ) since under time re- 
versal {x, V, C) (x, —V, — C) and therefore both J and Jfp 
change their signs [see arguments given after Eq. (fT3])1. Hence 
we finally get the same linear response result of Eq. ( fTSl l. 

The generalization to arbitrary dimensions is straightfor- 
ward and we outline the white noise case. We consider a 
d-dimensional hypercubic lattice with points represented by 
1 = (^1, ^2, ^d) where = 1, 2... with a ~ l,2,..,d. 
Let xi and vi be the d-dimensional displacement and velocity 
vectors respectively, of the particle at 1. Heat conduction is 
assumed to take place in the a = direction because of heat 
baths at temperature Tl and Tr that are attached to all lattice 
points on the two hypersurfaces 1^ = 1 and = N. The 
corresponding Langevin equations of motion are: 



mivi = fi + (5i„,i[j7i^ - Ti^vi] + (5;„,l[j7^^ - 



where 1 — {l^, 1'), so that 1' denotes points on a constant l^, 
hypersurface. The noise terms at different lattice points and in 
different directions are assumed to be uncorrelated, and satisfy 
the usual fluctuation-dissipation relations. 

Defining the layer energy e;^ = J^v the interlayer 

current ^ we find, following the same steps as in the 
ID case, the analogue of Eq. (fTZt with J replaced by J'^ ~ 



N- 



+1.1, 



and Jf, replaced by: 



ijlL-JNR)/^ = 



»7i'-v(i,i') 



(19) 



Writing the FP equation and acting with on the equilib- 
rium distribution gives: 



J: 
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*(Af,l') 



^From the forms of and JJ^, it is clear that we can repeat 
the arguments for the ID case which led to Eqs. ( ISIlll l. Hence 
we get Eq. ([15]) with J replaced by = J''/{N - 1). 

Fluid system coupled to Maxwell baths: We first consider 
a ID system of particles in a box of length L. The end parti- 
cles (1 and N) interact with baths at temperatures Tl and Tr 
respectively. Whenever the first particle hits the left wall it 
is reflected back with a random velocity chosen from the dis- 
tribution: Il{v) = mi Pl &{v) V exp[— /?LTOiti^/2] , with a 
similar rule at the right end. Otherwise the dynamics is Hamil- 
tonian. 

We find the FP current by noting that Jfp = 
{Af3)^^[dtP/ P]p=pg. There are two parts to the evolution 
of the phase space density: the Hamiltonian dynamics inside 
the system, and the effect of the heat baths. After a small time 
interval e, the phase space density P(x; v;t + e) is 

= /Jimie"^'^^'"!^! J^P{0, x' - v'e; -wq, v' - a'e; t)vodvo 

for xi < Vie 

= /3flmAre-5''«™""« /o°°P(x' - v'e, L; v' - a'e, vq; t)vodvo 
for XN > L + VNe 
= P(x — ve, V — ae, t) otherwise (20) 

where the primed variables in the first and second lines leave 
out particles 1 and N respectively. (Note that since < xi 
and a; AT < L, the conditions in the second and third lines 
imply wi > and ujv < 0.) 

If Tl = Tr ^ T, and P(x, v,<) = Po, the equilibrium 
phase space density for the temperature T, then the phase 
space density at time t + e is the same. Now if Tl.r — 
T ± AT/2, with P(x, v, t) stiU equal to Pq, then 



P(x;v;i + e) 



^0 



AT 
"2T 



{\l3mivl - l)e{vie 



{\(3mNv1i - l)6{xN - L - VNi) 



Po. 



(18) Dividing by e throughout and taking e ^ , we see that 



Jfp = --{^miVj^ - kBT)vi5{xi)9{vi) 

- -^UrnNvl, - kBT)vNS{xN - L)9{-vn)- (21) 

We have to use continuum energy density e(x, t) and current 
j{x,t), and the total heat current is now J = j j{x)dx in- 
stead of The continuity equation is still valid, and 

defining D{x,t) = J^dx'e{x',t) ~ J^e{x',t) and A{t) = 

dxD{x), we get the analogue of Eq. (flSIl : 



{J{t)J{0))dt^ L / {Jbit)J{0))dt. (22) 
Jo 
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Here Jf, = ^[ji.L — jN,R] as before, and 

ji.L = ^miviivfj^ - vf) 6{xi)9{-vi) 
3n,r = ^mNVN{v% L - vl) S{xN - L)0{vn) ■ 

The (5-functions enforce the condition that the particle is col- 
Hding with the bath, and vi,l and WAr.jj are the random ve- 
locities with which they emerge from the collision. Invoking 
detailed balance, using the explicit forms of J/p and Jh, and 
the fact that J(0) is uncorrelated with vi^lA'n.r we can show 
that {J{Q)Jb{t)) ^ -{J{t)Jfp{0)). With Eqs. ^ and (H), 
we obtain Eq. ( fTSl ) with {N — 1) replaced with L. The gener- 
alization to a (i-dimensional system is straightforward. First, 
any particle can interact with the baths at the ends if it reaches 
X — 01 X = L. Including the effect of the components of the 
velocity transverse to the heat-flow direction the derivation of 
Eq. ( I2T] ) gets modified and gives 

Jfp = Y.-l('2mvf - Ud+l)kBTKS{xr)9ivr) 

- - Ud+l)kBTXS{xr ~ m-v^). 

The expression for changes similarly, so that the final result 
of the previous paragraph is still valid. 

Conclusions: In this paper we have derived an exact expres- 
sion for the linear response conductance in a system connected 



to heat baths. Our results are valid in arbitrary dimensions 
and have been derived both for a solid where particles execute 
small displacements about fixed lattice positions as well for a 
fluid system where the motion of particles is unrestricted, and 
various heat bath models have been considered. 

The important differences with the usual Green-Kubo for- 
mula are worth noting. In the present formula, one does not 
need to first take the limit of infinite system size; the result 
is valid for finite systems. The fact that a sensible answer is 
obtained even for a finite system (unlike the case for the usual 
Green-Kubo formula) is because here we are dealing with an 
open system. Secondly the correlation function here has to 
be evaluated not with Hamiltonian dynamics, but for an open 
system evolving with heat bath dynamics. Finally we note 
that unlike the usual derivation of the Green-Kubo formula 
where the assumption of local thermal equilibrium is crucial, 
the present derivation requires no such assumption. The re- 
sults are thus valid even for integrable Hamiltonian models, 
the only requirement being that they should attain thermal 
equilibrium when coupled to one or more heat reservoirs all 
at the same temperature. 

Our derivation here is based on using both the microscopic 
equations of motion and also the equation for the phase space 
distribution. The broad class of systems and heat baths for 
which we have obtained our results strongly suggests that they 
are valid whenever detailed balance is satisfied. 
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